# Script produces figures presented in article and appendix  
# Last Updated: January 15, 2024

# tidyverse (Version 1.3.1)
if((!"tidyverse" %in% installed.packages()) | # If package is not installed 
   packageVersion("tidyverse") != "1.3.1"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/tidyverse/tidyverse_1.3.1.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("tidyverse")
} else { # Otherwise, load package
  library("tidyverse")
}

# stargazer (Version 5.2.1)
if((!"stargazer" %in% installed.packages()) | # If package is not installed 
   packageVersion("stargazer") != "5.2.1"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/stargazer/stargazer_5.2.1.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("stargazer")
} else { # Otherwise, load package
  library("stargazer")
}

# PanelMatch (Version 1.0.0)
if((!"PanelMatch" %in% installed.packages()) | # If package is not installed 
   packageVersion("PanelMatch") != "1.0.0"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/PanelMatch/PanelMatch_1.0.0.tar.gz" 
  install.packages(package_url, repos=NULL, type="source")
  library("PanelMatch")
} else { # Otherwise, load package
  library("PanelMatch")
}

# xtable (Version 1.8.3)
if((!"xtable" %in% installed.packages()) | # If package is not installed 
   packageVersion("xtable") != "1.8.3"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/xtable/xtable_1.8-3.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("xtable")
} else { # Otherwise, load package
  library("xtable")
}

# Load data and results
load("data/analysis_data.RData") # Analysis data
load("data/results.RData") # Results

# Figure 1 ----- 
graph_data <- data.frame()
for(i in 1:nrow(unilateral_exits)){
  dat_temp <- data_full |>
    filter(
      decision_year >= unilateral_exits$decision_year[i] - 10 &
        decision_year <= unilateral_exits$decision_year[i] + 10 &  
        ratification == 1 &  
        part_ccode == unilateral_exits$part_ccode[i]) |>
    mutate(
      t = decision_year - unilateral_exits$decision_year[i]
      ) |>
    group_by(t) |>
    summarize(
      y = sum(ratification, na.rm = T),
      n = 1, 
      state.b = unilateral_exits$part_ccode[i],
      w_treaty = unilateral_exits$treaty_id[i],
      w_year = unilateral_exits$decision_year[i]
      )
  graph_data <- rbind(graph_data, dat_temp)
}

# Average across cases
graph_data <- graph_data %>% group_by(t) %>%
  summarize(
    # Behavior of withdrawing state
    y_ci975 = mean(y, na.rm = T) + (qnorm(0.975) *(sd(y, na.rm = T)/sqrt(sum(n)))), # Get CIs
    y_ci025 = mean(y, na.rm = T) - (qnorm(0.975) *(sd(y, na.rm = T)/sqrt(sum(n)))),
    y = mean(y, na.rm = T) # Summarize Outcomes
    
  )

# Plot
ggplot(graph_data, aes(x = t, y = y)) +  
  scale_fill_brewer(palette= "Dark2") + 
  geom_line(aes(x = t, y = y), alpha = 1, linewidth = 0.75) + 
  geom_line(aes(x = t, y = y_ci975), color = "grey", linetype = "dashed", linewidth = 0.75) +
  geom_line(aes(x = t, y = y_ci025), color = "grey", linetype = "dashed", linewidth = 0.75) +
  theme_bw() +
  labs(x = "Year Relative to Exit (t = 0 is period of withdrawal)", 
       y = "Average Number of Multilateral Treaty \n Ratifications by Withdrawing States") +
  geom_vline(xintercept = seq(0, 15, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-10, 10)) + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=16),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=16), 
        text = element_text(size = 16),
        axis.title = element_text(size = 18), 
        axis.text = element_text(size = 14)) 
ggsave("figs/isolationism.jpeg", plot = last_plot(), width = 8, height = 5, units = "in")

# Import regime ideal point values
load("data/preprocessing/AgreementScoresAll_Jul2023.Rdata") # Import years of UNGA sessions
dfAgree <- dfAgree |> group_by(session.x) |> summarize(year = mean(year, na.rm = T))
colnames(dfAgree) <- c("session", "year")
ideal_points <- read.csv("data/preprocessing/IdealpointestimatesAll_Jul2023.csv") # Import UNGA Ideal Point data 
ideal_points <- left_join(ideal_points, dfAgree, by = "session") # Merge

graph_data <- data.frame()
for(i in 1:nrow(unilateral_exits)){
  dat_temp <- ideal_points |>
    filter(
      year >= unilateral_exits$decision_year[i] - 10 &
        year <= unilateral_exits$decision_year[i] + 10 &  
        ccode == unilateral_exits$part_ccode[i]) |>
    mutate(
      t = year - unilateral_exits$decision_year[i]
    ) |>
    group_by(t) |>
    summarize(
      y = mean(IdealPointAll, na.rm = T),
      n = 1, 
      state.b = unilateral_exits$part_ccode[i],
      w_treaty = unilateral_exits$treaty_id[i],
      w_year = unilateral_exits$decision_year[i]
    )
  graph_data <- rbind(graph_data, dat_temp)
}

graph_data <- graph_data %>% group_by(t) %>%
  summarize(
    # Behavior of withdrawing state
    y_ci975 = mean(y, na.rm = T) + (qnorm(0.975) *(sd(y, na.rm = T)/sqrt(sum(n)))), # Get CIs
    y_ci025 = mean(y, na.rm = T) - (qnorm(0.975) *(sd(y, na.rm = T)/sqrt(sum(n)))),
    y = mean(y, na.rm = T) # Summarize Outcomes
    
  )

# Plot
ggplot(graph_data, aes(x = t, y = y)) +  
  scale_fill_brewer(palette= "Dark2") + 
  geom_line(aes(x = t, y = y), alpha = 1, linewidth = 0.75) + 
  geom_line(aes(x = t, y = y_ci975), color = "grey", linetype = "dashed", linewidth = 0.75) +
  geom_line(aes(x = t, y = y_ci025), color = "grey", linetype = "dashed", linewidth = 0.75) +
  theme_bw() +
  labs(x = "Year Relative to Exit (t = 0 is period of withdrawal)", 
       y = "Average UNGA Ideal Point \n of Withdrawing States") +
  geom_vline(xintercept = seq(0, 15, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-10, 10)) + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=16),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=16), 
        text = element_text(size = 16),
        axis.title = element_text(size = 18), 
        axis.text = element_text(size = 14)) 
ggsave("figs/unga_preference_change.jpeg", plot = last_plot(), width = 8, height = 5, units = "in")

# Create function to extract information from PanelMatch objects
# Set confidence levels
high_ci <- 0.05 
low_ci <- 0.1

# Make simplifying functions
high_ci_upper_bound_fxn <- function(x){quantile(x, probs = c(1-high_ci/2))}
high_ci_lower_bound_fxn <- function(x){quantile(x, probs = c(high_ci/2))}
low_ci_upper_bound_fxn <- function(x){quantile(x, probs = c(1-low_ci/2))}
low_ci_lower_bound_fxn <- function(x){quantile(x, probs = c(low_ci/2))}

pm_extractor_fxn <- function(analysis, sample, ests, upper_critical_value, lower_critical_value){
  est_names <- names(ests$estimates)
  est_betas <- unname(ests$estimates)
  est_ses <- ests$standard.error
  boots <- ests$bootstrapped.estimates
  
  dat_temp <- c(rep(analysis, length(est_names)), 
                rep(sample, length(est_names)),
                est_betas,
                est_ses,
                apply(boots, 2, low_ci_lower_bound_fxn),
                apply(boots, 2, low_ci_upper_bound_fxn),
                apply(boots, 2, high_ci_lower_bound_fxn),
                apply(boots, 2, high_ci_upper_bound_fxn),
                est_names)
  
  table_out <- matrix(data = dat_temp, nrow = length(est_ses), ncol = 9, byrow = F)
  colnames(table_out) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
  return(table_out)
}

# Direct Effects ----
# Extract true estimates and confidence intervals
graph_values <- as.data.frame(matrix(NA, ncol = 9, nrow = 0))
colnames(graph_values) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "All Exits", ests_cbps_direct, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Major Power Exits", ests_cbps_direct_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Minor Power Exits", ests_cbps_direct_not_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 1st Tertile", ests_cbps_direct_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 2nd Tertile", ests_cbps_direct_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 3rd Tertile", ests_cbps_direct_tertile3, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "All Exits (Restricted)", ests_cbps_direct_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Major Power Exits (Restricted)", ests_cbps_direct_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Minor Power Exits (Restricted)", ests_cbps_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 1st Tertile (Restricted)", ests_cbps_direct_tertile1_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 2nd Tertile (Restricted)", ests_cbps_direct_tertile2_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 3rd Tertile (Restricted)", ests_cbps_direct_tertile3_restricted, high_ci, low_ci))

graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "All Exits", ests_maha5_direct, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Major Power Exits", ests_maha5_direct_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Minor Power Exits", ests_maha5_direct_not_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 1st Tertile", ests_maha5_direct_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 2nd Tertile", ests_maha5_direct_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 3rd Tertile", ests_maha5_direct_tertile3, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "All Exits (Restricted)", ests_maha5_direct_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Major Power Exits (Restricted)", ests_maha5_direct_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Minor Power Exits (Restricted)", ests_maha5_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 1st Tertile (Restricted)", ests_maha5_direct_tertile1_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 2nd Tertile (Restricted)", ests_maha5_direct_tertile2_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 3rd Tertile (Restricted)", ests_maha5_direct_tertile3_restricted, high_ci, low_ci))

graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "All Exits", ests_noref_direct, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Major Power Exits", ests_noref_direct_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Minor Power Exits", ests_noref_direct_not_rel_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 1st Tertile", ests_noref_direct_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 2nd Tertile", ests_noref_direct_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 3rd Tertile", ests_noref_direct_tertile3, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "All Exits (Restricted)", ests_noref_direct_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Major Power Exits (Restricted)", ests_noref_direct_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Minor Power Exits (Restricted)", ests_noref_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 1st Tertile (Restricted)", ests_noref_direct_tertile1_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 2nd Tertile (Restricted)", ests_noref_direct_tertile2_restricted, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 3rd Tertile (Restricted)", ests_noref_direct_tertile3_restricted, high_ci, low_ci))

# Extract placebo estimates and confidence intervals
graph_values_placebos <- as.data.frame(matrix(NA, ncol = 9, nrow = 0))
colnames(graph_values_placebos) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "All Exits", ests_cbps_p2_direct, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Major Power Exits", ests_cbps_p2_direct_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Minor Power Exits", ests_cbps_p2_direct_not_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 1st Tertile", ests_cbps_p2_direct_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 2nd Tertile", ests_cbps_p2_direct_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 3rd Tertile", ests_cbps_p2_direct_tertile3, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "All Exits (Restricted)", ests_cbps_p2_direct_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Major Power Exits (Restricted)", ests_cbps_p2_direct_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Minor Power Exits (Restricted)", ests_cbps_p2_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 1st Tertile (Restricted)", ests_cbps_p2_direct_tertile1_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 2nd Tertile (Restricted)", ests_cbps_p2_direct_tertile2_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 3rd Tertile (Restricted)", ests_cbps_p2_direct_tertile3_restricted, high_ci, low_ci))

graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "All Exits", ests_maha5_p2_direct, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Major Power Exits", ests_maha5_p2_direct_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Minor Power Exits", ests_maha5_p2_direct_not_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 1st Tertile", ests_maha5_p2_direct_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 2nd Tertile", ests_maha5_p2_direct_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 3rd Tertile", ests_maha5_p2_direct_tertile3, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "All Exits (Restricted)", ests_maha5_p2_direct_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Major Power Exits (Restricted)", ests_maha5_p2_direct_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Minor Power Exits (Restricted)", ests_maha5_p2_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 1st Tertile (Restricted)", ests_maha5_p2_direct_tertile1_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 2nd Tertile (Restricted)", ests_maha5_p2_direct_tertile2_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 3rd Tertile (Restricted)", ests_maha5_p2_direct_tertile3_restricted, high_ci, low_ci))

graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "All Exits", ests_noref_p2_direct, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Major Power Exits", ests_noref_p2_direct_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Minor Power Exits", ests_noref_p2_direct_not_rel_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 1st Tertile", ests_noref_p2_direct_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 2nd Tertile", ests_noref_p2_direct_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 3rd Tertile", ests_noref_p2_direct_tertile3, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "All Exits (Restricted)", ests_noref_p2_direct_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Major Power Exits (Restricted)", ests_noref_p2_direct_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Minor Power Exits (Restricted)", ests_noref_p2_direct_not_rel_power_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 1st Tertile (Restricted)", ests_noref_p2_direct_tertile1_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 2nd Tertile (Restricted)", ests_noref_p2_direct_tertile2_restricted, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 3rd Tertile (Restricted)", ests_noref_p2_direct_tertile3_restricted, high_ci, low_ci))

# Fix time series 
graph_values_placebos$t[graph_values_placebos$t == "t+0"] <- -2
graph_values_placebos$t[graph_values_placebos$t == "t+1"] <- -1
graph_values$t[graph_values$t == "t+0"] <- 0
graph_values$t[graph_values$t == "t+1"] <- 1
graph_values$t[graph_values$t == "t+2"] <- 2
graph_values$t[graph_values$t == "t+3"] <- 3
graph_values$t[graph_values$t == "t+4"] <- 4
graph_values$t[graph_values$t == "t+5"] <- 5
graph_values$t[graph_values$t == "t+6"] <- 6

# Combine true and placebo estimates
graph_values <- rbind(graph_values, graph_values_placebos)

for(i in 3:9){ # Coerce estimates and SEs to numeric
  graph_values[,i] <- as.numeric(graph_values[,i])
}

# Make function to automate plots
plot_fxn <- function(x,y){
  out <- graph_values %>%
    filter(pooling == x & name == y) %>%
    ggplot(aes(x = t, y = estimate)) +  
    scale_color_brewer(palette= "Dark2") + 
    geom_crossbar(aes(ymin = lower_95, ymax = upper_95), 
                  width=0.1, fill = "black", color = "black") +
    geom_point(aes(x = t, y = estimate), size = 5, alpha = 1) + 
    scale_shape_manual(values=c(16, 1)) + 
    theme_bw() +
    geom_hline(yintercept=0, linetype="dashed", color = "black") +
    geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
    coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.75, 0.75)) + 
    scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                       labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
    geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
    ylab("Change in # reforms") +
    xlab("Years relative to unilateral exit") + 
    theme(legend.position= "none", 
          legend.title=element_blank(), 
          legend.text = element_text(size=14),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 14),
          plot.title = element_text(hjust = 0.5, size=16), 
          text = element_text(size = 20),
          axis.title = element_text(size = 20), 
          axis.text = element_text(size = 16)) 
  return(out)
}

# Supp. Coefficient Plots of Direct Effect -----
plot_fxn(y = "All Exits", x = "CBPS")
ggsave("figs/direct_cbps.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "CBPS")
ggsave("figs/direct_cbps_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "CBPS")
ggsave("figs/direct_cbps_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits", x = "Maha")
ggsave("figs/direct_maha.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "Maha")
ggsave("figs/direct_maha_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "Maha")
ggsave("figs/direct_maha_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits", x = "NoRef")
ggsave("figs/direct_noref.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "NoRef")
ggsave("figs/direct_noref_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "NoRef")
ggsave("figs/direct_noref_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits (Restricted)", x = "Maha")
ggsave("figs/direct_maha_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits (Restricted)", x = "NoRef")
ggsave("figs/direct_noref_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits (Restricted)", x = "Maha")
ggsave("figs/direct_maha_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits (Restricted)", x = "NoRef")
ggsave("figs/direct_noref_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_not_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits (Restricted)", x = "Maha")
ggsave("figs/direct_maha_not_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits (Restricted)", x = "NoRef")
ggsave("figs/direct_noref_not_power_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 1st Tertile", x = "CBPS")
ggsave("figs/direct_cbps_tertile1.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 2nd Tertile", x = "CBPS")
ggsave("figs/direct_cbps_tertile2.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 3rd Tertile", x = "CBPS")
ggsave("figs/direct_cbps_tertile3.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 1st Tertile (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_tertile1_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 2nd Tertile (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_tertile2_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 3rd Tertile (Restricted)", x = "CBPS")
ggsave("figs/direct_cbps_tertile3_restricted.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

# Figure 4 ----
direct_main_out <- graph_values %>% 
  filter(pooling == "CBPS" & name == "All Exits") %>% 
  ggplot(aes(x = t, y = estimate)) +  
  geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
              width=0.075, position = position_dodge(width = 0.75)) +
  geom_crossbar(aes(ymin = lower_90, ymax = upper_90), show.legend = F,
                width=0.05, position = position_dodge(width = 0.75)) +
  geom_point(aes(x = t, y = estimate),
  size = 3, show.legend = F,
  position = position_dodge(width = 0.75)) +
  scale_shape_manual(values=c(1, 16)) + 
  scale_color_brewer(palette= "Dark2") + 
  theme_bw() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.25, 0.5)) + 
  scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                     labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
  geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
  ylab("Change in number of reforms (ATT)") +
  xlab("Years relative to unilateral exit (t)") + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=12),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=14), 
        text = element_text(size = 12),
        axis.title = element_text(size = 12), 
        axis.text = element_text(size = 12)) 
direct_main_out
ggsave("figs/direct_main.jpeg", plot = last_plot(), width = 5, height = 3.25, units = "in")

# Figure 5 ----
direct_main_power_out <- graph_values %>% 
  filter(pooling == "CBPS" & 
           (name == "Major Power Exits" | 
              name == "Minor Power Exits")) %>% 
  ggplot(aes(x = t, y = estimate)) +  
  facet_wrap(~ name, nrow=1, scales = "fixed",  strip.position = c("top")) + 
  geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
              width=0.075, position = position_dodge(width = 0.75)) +
  geom_crossbar(aes(ymin = lower_90, ymax = upper_90), show.legend = F,
                width=0.05, position = position_dodge(width = 0.75)) +
  geom_point(aes(x = t, y = estimate),
  size = 3, show.legend = F,
  position = position_dodge(width = 0.75)) +
  scale_shape_manual(values=c(1, 16)) + 
  theme_bw() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.5, 0.75)) + 
  scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                     labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
  geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
  ylab("Change in number of reforms (ATT)") +
  xlab("Years relative to unilateral exit (t)") + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=14),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=14), 
        text = element_text(size = 14),
        axis.title = element_text(size = 14), 
        axis.text = element_text(size = 14)) 
direct_main_power_out
ggsave("figs/direct_power_main.jpeg", plot = last_plot(), width = 9, height = 4, units = "in")

# Indirect Effects ----
graph_values <- as.data.frame(matrix(NA, ncol = 9, nrow = 0))
colnames(graph_values) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "All Exits", ests_cbps_indirect, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Major Power Exits", ests_cbps_indirect_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Minor Power Exits", ests_cbps_indirect_not_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 1st Tertile", ests_cbps_indirect_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 2nd Tertile", ests_cbps_indirect_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("CBPS", "Power 3rd Tertile", ests_cbps_indirect_tertile3, high_ci, low_ci))

graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "All Exits", ests_maha5_indirect, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Major Power Exits", ests_maha5_indirect_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Minor Power Exits", ests_maha5_indirect_not_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 1st Tertile", ests_maha5_indirect_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 2nd Tertile", ests_maha5_indirect_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("Maha", "Power 3rd Tertile", ests_maha5_indirect_tertile3, high_ci, low_ci))

graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "All Exits", ests_noref_indirect, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Major Power Exits", ests_noref_indirect_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Minor Power Exits", ests_noref_indirect_not_power, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 1st Tertile", ests_noref_indirect_tertile1, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 2nd Tertile", ests_noref_indirect_tertile2, high_ci, low_ci))
graph_values <- rbind(graph_values, pm_extractor_fxn("NoRef", "Power 3rd Tertile", ests_noref_indirect_tertile3, high_ci, low_ci))

# Combine placebos
graph_values_placebos <- as.data.frame(matrix(NA, ncol = 9, nrow = 0))
colnames(graph_values_placebos) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "All Exits", ests_cbps_p2_indirect, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Major Power Exits", ests_cbps_p2_indirect_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Minor Power Exits", ests_cbps_p2_indirect_not_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 1st Tertile", ests_cbps_p2_indirect_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 2nd Tertile", ests_cbps_p2_indirect_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("CBPS", "Power 3rd Tertile", ests_cbps_p2_indirect_tertile3, high_ci, low_ci))

graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "All Exits", ests_maha5_p2_indirect, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Major Power Exits", ests_maha5_p2_indirect_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Minor Power Exits", ests_maha5_p2_indirect_not_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 1st Tertile", ests_maha5_p2_indirect_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 2nd Tertile", ests_maha5_p2_indirect_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("Maha", "Power 3rd Tertile", ests_maha5_p2_indirect_tertile3, high_ci, low_ci))

graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "All Exits", ests_noref_p2_indirect, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Major Power Exits", ests_noref_p2_indirect_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Minor Power Exits", ests_noref_p2_indirect_not_power, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 1st Tertile", ests_noref_p2_indirect_tertile1, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 2nd Tertile", ests_noref_p2_indirect_tertile2, high_ci, low_ci))
graph_values_placebos <- rbind(graph_values_placebos, pm_extractor_fxn("NoRef", "Power 3rd Tertile", ests_noref_p2_indirect_tertile3, high_ci, low_ci))

# Fix time series 
graph_values_placebos$t[graph_values_placebos$t == "t+0"] <- -2
graph_values_placebos$t[graph_values_placebos$t == "t+1"] <- -1
graph_values$t[graph_values$t == "t+0"] <- 0
graph_values$t[graph_values$t == "t+1"] <- 1
graph_values$t[graph_values$t == "t+2"] <- 2
graph_values$t[graph_values$t == "t+3"] <- 3
graph_values$t[graph_values$t == "t+4"] <- 4
graph_values$t[graph_values$t == "t+5"] <- 5
graph_values$t[graph_values$t == "t+6"] <- 6

# Combine true and placebo estimates
graph_values <- rbind(graph_values, graph_values_placebos)

for(i in 3:9){ # Coerce estimates and SEs to numeric
  graph_values[,i] <- as.numeric(graph_values[,i])
}

# Make function to automate plots
plot_fxn <- function(x,y){
  out <- graph_values %>%
    filter(pooling == x & name == y) %>%
    ggplot(aes(x = t, y = estimate)) +  
    scale_color_brewer(palette= "Dark2") + 
    geom_crossbar(aes(ymin = lower_95, ymax = upper_95), 
                  width=0.1, fill = "black", color = "black") +
    geom_point(aes(x = t, y = estimate), size = 5) + 
    scale_shape_manual(values=c(16, 1)) + 
    theme_bw() +
    geom_hline(yintercept=0, linetype="dashed", color = "black") +
    geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
    coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.125, 0.125) ) + #
    scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                       labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
    geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
    ylab("Change in # reforms") +
    xlab("Years relative to unilateral exit") + 
    theme(legend.position= "none", 
          legend.title=element_blank(), 
          legend.text = element_text(size=14),
          strip.text = element_text(size = 14),
          panel.grid.minor = element_blank(),
          plot.title = element_text(hjust = 0.5, size=16), 
          text = element_text(size = 20),
          axis.title = element_text(size = 20), 
          axis.text = element_text(size = 16)) 
  return(out)
}

# Supp. Coefficient Plots of Indirect Effect -----
plot_fxn(y = "All Exits", x = "CBPS")
ggsave("figs/indirect_cbps.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "CBPS")
ggsave("figs/indirect_cbps_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "CBPS")
ggsave("figs/indirect_cbps_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits", x = "Maha")
ggsave("figs/indirect_maha.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "Maha")
ggsave("figs/indirect_maha_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "Maha")
ggsave("figs/indirect_maha_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "All Exits", x = "NoRef")
ggsave("figs/indirect_noref.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Major Power Exits", x = "NoRef")
ggsave("figs/indirect_noref_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Minor Power Exits", x = "NoRef")
ggsave("figs/indirect_noref_not_power.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 1st Tertile", x = "CBPS")
ggsave("figs/indirect_cbps_tertile1.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 2nd Tertile", x = "CBPS")
ggsave("figs/indirect_cbps_tertile2.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")

plot_fxn(y = "Power 3rd Tertile", x = "CBPS")
ggsave("figs/indirect_cbps_tertile3.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")


# Figure 6 ----
indirect_main_out <- graph_values %>% 
  filter(pooling == "CBPS" & name == "All Exits") %>% 
  ggplot(aes(x = t, y = estimate)) +  
  geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
              width=0.075, position = position_dodge(width = 0.75)) +
  geom_crossbar(aes(ymin = lower_90, ymax = upper_90), show.legend = F,
                width=0.05, position = position_dodge(width = 0.75)) +
  geom_point(aes(x = t, y = estimate),
  size = 3, show.legend = F,
  position = position_dodge(width = 0.75)) +
  scale_shape_manual(values=c(1, 16)) + 
  scale_color_brewer(palette= "Dark2") + 
  theme_bw() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.1, 0.05)) + 
  scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                     labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
  geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
  ylab("Change in number of reforms (ATT)") +  
  xlab("Years relative to unilateral exit (t)") + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=12),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=14), 
        text = element_text(size = 12),
        axis.title = element_text(size = 12), 
        axis.text = element_text(size = 12)) 
indirect_main_out
ggsave("figs/indirect_main.jpeg", plot = last_plot(), width = 5, height = 3.25, units = "in")

# Figure 7 ----
indirect_main_power_out <- 
  graph_values |> 
  filter(pooling == "CBPS" & 
        (name == "Major Power Exits" | 
         name == "Minor Power Exits")) |> 
  ggplot(aes(x = t, y = estimate)) +  
  facet_wrap(~ name, nrow=1, scales = "fixed",  strip.position = c("top")) + 
  geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
              width=0.075, position = position_dodge(width = 0.75)) +
  geom_crossbar(aes(ymin = lower_90, ymax = upper_90), show.legend = F,
                width=0.05, position = position_dodge(width = 0.75)) +
  geom_point(aes(x = t, y = estimate),
  size = 3, show.legend = F,
  position = position_dodge(width = 0.75)) +
  scale_shape_manual(values=c(1, 16)) + 
  theme_bw() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept = seq(-0.5, 10, by = 0.01), colour="grey", linetype = "solid", alpha = 0.05) +
  coord_cartesian(xlim = c(-2.25, 4.25), ylim = c(-0.125, 0.075)) + 
  scale_x_continuous(breaks=c(-5, -4, -3, -2, -1, 0,1, 2, 3, 4, 5, 6),
                     labels=c("-5", "-4", "-3", "-2", "-1", "0", "1", "2", "3", "4", "5", "6")) + 
  geom_rect(aes(xmin = -0.5, xmax = 10, ymin = -1, ymax = 1), colour="grey", alpha = 0.02) + 
  ylab("Change in number of reforms (ATT)") + 
  xlab("Years relative to unilateral exit (t)") + 
  theme(legend.position= "bottom", 
        legend.title=element_blank(), 
        legend.text = element_text(size=14),
        strip.text = element_text(size = 16),
        plot.title = element_text(hjust = 0.5, size=14), 
        text = element_text(size = 14),
        axis.title = element_text(size = 14), 
        axis.text = element_text(size = 14)) 
indirect_main_power_out
ggsave("figs/indirect_power_main.jpeg", plot = last_plot(), width = 9, height = 4, units = "in")

# Balance Plots ----

# Function to automate balance plots
balance_plot_fxn <- function(x){
  # Extract standardized differences
  # Combine and format dataframe for plotting
  m_out <- as.data.frame(x)
  m_out$t <- -3:0
  m_out$DV <- x[, "reform_n"]
  m_out$analysis <- "Logged" # Names are just used a placeholders
  m_out$refinement <- "Mahalanobis Distance \n (5 Matches)" 
  graph_values_long <- m_out %>% pivot_longer( # Make long format
    cols = all_of(cols_to_make_long),
    values_to = "values")
  graph_values_long$refinement <- fct_rev(graph_values_long$refinement)
  graph_values_long <- subset(graph_values_long, t < 0) # Drop t = 0, which is post-treatment
  graph_values_long_dv <- graph_values_long %>% filter(name == "reform_n")
  graph_values_long$refinement <- factor(graph_values_long$refinement) # Reorder factors 
  plot_out <- ggplot(graph_values_long, # Make Plot
                     aes(x = t, y = values, group = name, linetype = 1)) +  
    scale_fill_brewer(palette= "Dark2") + 
    geom_line(aes(x = t, y = values, linetype = name, group = name), alpha = 0.25, linewidth = 0.75) + 
    geom_line(data = graph_values_long_dv, col = "black", alpha = 1, linetype = 1, linewidth = 1) +
    theme_bw() +
    scale_x_continuous(breaks=c(-3, -2, -1),
                       labels=c( "-3", "-2", "-1")) + 
    labs(x = "Years relative to exit", 
         y = "Standardized Mean Diff.") +
    geom_hline(yintercept = 0, colour="black", linetype = "dashed", alpha = 1) +
    coord_cartesian(ylim = c(-1, 1)) +
    theme(legend.position= "none", 
          legend.title=element_blank(), 
          strip.text = element_text(size = 10),
          legend.text = element_text(size=14),
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank(),
          panel.grid.minor.x = element_blank(),
          plot.title = element_text(hjust = 0.5, size=16), 
          text = element_text(size = 20),
          axis.title = element_text(size = 20), 
          axis.text = element_text(size = 20))
  return(plot_out)
}

# Supp. Balance Plots -----
# List of variables used in all balance plots
balance_var_list <- c("reform_n", "n_members_log", 
                      "last_reform", "cumulative_agreements",
                      "pref_mean", "pref_sd",
                      "human_rights", "environment","security", "economics") 

# Indirect Balance Plots
cols_to_make_long  <- c("overlap", balance_var_list) # Assign relevant overlap variable to var list
balance_plot_fxn(balance_noref_indirect) # Make plot
ggsave("figs/balance_plot_indirect_noref.jpeg", plot = last_plot(), width = 4, height = 4, units = "in") # Save
balance_plot_fxn(balance_maha5_indirect) # Repeat
ggsave("figs/balance_plot_indirect_maha5.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_indirect)
ggsave("figs/balance_plot_indirect_cbps.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Indirect Power Balance Plots
cols_to_make_long  <- c("overlap_power", balance_var_list) 
balance_plot_fxn(balance_noref_indirect_power)
ggsave("figs/balance_plot_indirect_noref_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_indirect_power)
ggsave("figs/balance_plot_indirect_maha5_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_indirect_power)
ggsave("figs/balance_plot_indirect_cbps_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Indirect Not Power Balance Plots
cols_to_make_long  <- c("overlap_not_power", balance_var_list) 
balance_plot_fxn(balance_noref_indirect_not_power)
ggsave("figs/balance_plot_indirect_noref_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_indirect_not_power)
ggsave("figs/balance_plot_indirect_maha5_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_indirect_not_power)
ggsave("figs/balance_plot_indirect_cbps_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Balance Plots 
cols_to_make_long  <- c("overlap", balance_var_list) 
balance_plot_fxn(balance_noref_direct)
ggsave("figs/balance_plot_direct_noref.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct)
ggsave("figs/balance_plot_direct_maha5.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct)
ggsave("figs/balance_plot_direct_cbps.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct (Restricted) Balance Plots 
cols_to_make_long  <- c("overlap", balance_var_list) 
balance_plot_fxn(balance_noref_direct_restricted)
ggsave("figs/balance_plot_direct_noref_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct_restricted)
ggsave("figs/balance_plot_direct_maha5_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_restricted)
ggsave("figs/balance_plot_direct_cbps_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Power (Restricted) Balance Plots  
cols_to_make_long  <- c("overlap_rel_power", balance_var_list) 
balance_plot_fxn(balance_noref_direct_rel_power_restricted)
ggsave("figs/balance_plot_direct_noref_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct_rel_power_restricted)
ggsave("figs/balance_plot_direct_maha5_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_rel_power_restricted)
ggsave("figs/balance_plot_direct_cbps_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Not Power (Restricted) Balance Plots  
cols_to_make_long  <- c("overlap_not_rel_power", balance_var_list) 
balance_plot_fxn(balance_noref_direct_not_rel_power_restricted)
ggsave("figs/balance_plot_direct_noref_not_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct_not_rel_power_restricted)
ggsave("figs/balance_plot_direct_maha5_not_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_not_rel_power_restricted)
ggsave("figs/balance_plot_direct_cbps_not_power_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Power Balance Plots  
cols_to_make_long  <- c("overlap_rel_power", balance_var_list) 
balance_plot_fxn(balance_noref_direct_rel_power)
ggsave("figs/balance_plot_direct_noref_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct_rel_power)
ggsave("figs/balance_plot_direct_maha5_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_rel_power)
ggsave("figs/balance_plot_direct_cbps_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Not Power (Restricted) Balance Plots  
cols_to_make_long  <- c("overlap_not_rel_power", balance_var_list)
balance_plot_fxn(balance_noref_direct_not_rel_power)
ggsave("figs/balance_plot_direct_noref_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_maha5_direct_not_rel_power)
ggsave("figs/balance_plot_direct_maha5_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_not_rel_power)
ggsave("figs/balance_plot_direct_cbps_not_power.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Power Tertile 1 Balance Plots 
cols_to_make_long  <- c("overlap_rel_power_tertile1", balance_var_list)
balance_plot_fxn(balance_cbps_direct_tertile1)
ggsave("figs/balance_plot_direct_tertile1.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_tertile1_restricted)
ggsave("figs/balance_plot_direct_tertile1_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Indirect Power Tertile 1 Balance Plots 
cols_to_make_long  <- c("overlap_power_tertile1", balance_var_list)
balance_plot_fxn(balance_cbps_indirect_tertile1)
ggsave("figs/balance_plot_indirect_tertile1.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Power Tertile 2 Balance Plots 
cols_to_make_long  <- c("overlap_rel_power_tertile2", balance_var_list)
balance_plot_fxn(balance_cbps_direct_tertile2)
ggsave("figs/balance_plot_direct_tertile2.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_tertile2_restricted)
ggsave("figs/balance_plot_direct_tertile2_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Indirect Power Tertile 2 Balance Plots 
cols_to_make_long  <- c("overlap_power_tertile2", balance_var_list)
balance_plot_fxn(balance_cbps_indirect_tertile2)
ggsave("figs/balance_plot_indirect_tertile2.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Direct Power Tertile 3 Balance Plots 
cols_to_make_long  <- c("overlap_rel_power_tertile3", balance_var_list)
balance_plot_fxn(balance_cbps_direct_tertile3)
ggsave("figs/balance_plot_direct_tertile3.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")
balance_plot_fxn(balance_cbps_direct_tertile3_restricted)
ggsave("figs/balance_plot_direct_tertile3_restricted.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Indirect Power Tertile 3 Balance Plots 
cols_to_make_long  <- c("overlap_power_tertile3", balance_var_list)
balance_plot_fxn(balance_cbps_indirect_tertile3)
ggsave("figs/balance_plot_indirect_tertile3.jpeg", plot = last_plot(), width = 4, height = 4, units = "in")

# Figure 3 ----
m_out <- as.data.frame(balance_cbps_direct) # Extract balance data for main analyses
m_out$t <- -3:0
m_out$DV <- balance_cbps_direct[, "reform_n"]
m_out$analysis <- "Direct Effect"
m_out$refinement <- "CBPS Weighting"

m_out4 <- as.data.frame(balance_noref_direct)
m_out4$t <- -3:0
m_out4$DV <- balance_noref_direct[, "reform_n"]
m_out4$analysis <- "Direct Effect"
m_out4$refinement <- "No Refinement"

m_out2 <- as.data.frame(balance_cbps_indirect)
m_out2$t <- -3:0
m_out2$DV <- balance_cbps_indirect[, "reform_n"]
m_out2$analysis <- "Spillover Effect"
m_out2$refinement <- "CBPS Weighting"

m_out3 <- as.data.frame(balance_noref_indirect)
m_out3$t <- -3:0
m_out3$DV <- balance_cbps_indirect[, "reform_n"]
m_out3$analysis <- "Spillover Effect"
m_out3$refinement <- "No Refinement"

cols_to_make_long  <- c("overlap", balance_var_list) # Variable names
m_out <- rbind(m_out, m_out2, m_out3, m_out4) # Combine data
graph_values_long <- m_out %>% pivot_longer( # Make long format
  cols = all_of(cols_to_make_long),
  values_to = "values")
graph_values_long$analysis <- fct_rev(graph_values_long$analysis) # Coerce to factor
graph_values_long <- subset(graph_values_long, t < 0) # Drop t = 0, which is post-treatment
graph_values_long$analysis <- factor(graph_values_long$analysis, # Order factor variables
                                     levels = c("Direct Effect", "Spillover Effect"))
graph_values_long$refinement <- factor(graph_values_long$refinement, 
                                       labels = c("No Refinement", "CBPS Weighting"),
                                       levels = c("No Refinement", "CBPS Weighting"))
graph_values_long_dv <- graph_values_long %>% filter(name == "reform_n") # Extract outcome values
plot_out <- ggplot(graph_values_long, # Make Plot
                   aes(x = t, y = values, group = name, linetype = 1)) +  
  facet_grid(vars(analysis), vars(refinement)) + 
  scale_fill_brewer(palette= "Dark2") + 
  geom_line(aes(x = t, y = values, linetype = name, group = name), alpha = 0.25, linewidth = 0.75) + 
  geom_line(data = graph_values_long_dv, col = "black", alpha = 1, linetype = 1, linewidth = 1) +
  theme_bw() +
  scale_x_continuous(breaks=c(-3, -2, -1),
                     labels=c( "-3", "-2", "-1")) + 
  labs(x = "Years relative to exit", 
       y = "Standardized Mean Difference") +
  geom_hline(yintercept = 0, colour="black", linetype = "dashed", alpha = 1) +
  coord_cartesian(ylim = c(-1, 1)) +
  theme(legend.position= "none", 
        legend.title=element_blank(), 
        strip.text = element_text(size = 14),
        legend.text = element_text(size=14),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        plot.title = element_text(hjust = 0.5, size=16), 
        text = element_text(size = 14),
        axis.title = element_text(size = 14), 
        axis.text = element_text(size = 14))
plot_out
ggsave("figs/balance_plot_main.jpeg", plot = last_plot(), width = 5, height = 4, units = "in")


# Covariate balance reported in manuscript
# Direct effect analysis
max(abs(balance_cbps_direct[1:3,"reform_n"]))
max(abs(balance_cbps_direct[1:3,]))

# Indirect effect analysis
max(abs(balance_cbps_indirect[1:3,"reform_n"]))
max(abs(balance_cbps_indirect[1:3,]))

# Placebos ----

# Function to automate placebo plots
placebo_plot_fxn <- function(p1,p2,p3,p4,p5){
  graph_values <- as.data.frame(matrix(NA, ncol = 9, nrow = 0))
  colnames(graph_values) <- c("pooling", "name", "estimate", "standard_error", "lower_90", "upper_90", "lower_95", "upper_95", "t")
  graph_values <- rbind(graph_values, pm_extractor_fxn("Logged", "t-1", p1, high_ci, low_ci))
  graph_values <- rbind(graph_values, pm_extractor_fxn("Logged", "t-2", p2, high_ci, low_ci))
  graph_values <- rbind(graph_values, pm_extractor_fxn("Logged", "t-3", p3, high_ci, low_ci))
  graph_values <- rbind(graph_values, pm_extractor_fxn("Logged", "t-4", p4, high_ci, low_ci))
  graph_values <- rbind(graph_values, pm_extractor_fxn("Logged", "t-5", p5, high_ci, low_ci))
  
  for(i in 3:8){ # Coerce estimates and SEs to numeric
    graph_values[,i] <- as.numeric(graph_values[,i])
  }
  
  # Correct t so that it is relative to start date instead of placebo
  graph_values$t[graph_values$name == "t-1"] <- "t+4" # Fix t-1
  graph_values$t[graph_values$name == "t-2" & graph_values$t == "t+1"] <- "t+4" # Fix t-2
  graph_values$t[graph_values$name == "t-2" & graph_values$t == "t+0"] <- "t+3" # Fix t-2
  graph_values$t[graph_values$name == "t-3" & graph_values$t == "t+2"] <- "t+4" # Fix t-3
  graph_values$t[graph_values$name == "t-3" & graph_values$t == "t+1"] <- "t+3" # Fix t-3
  graph_values$t[graph_values$name == "t-3" & graph_values$t == "t+0"] <- "t+2" # Fix t-3
  graph_values$t[graph_values$name == "t-4" & graph_values$t == "t+3"] <- "t+4" # Fix t-4
  graph_values$t[graph_values$name == "t-4" & graph_values$t == "t+2"] <- "t+3" # Fix t-4
  graph_values$t[graph_values$name == "t-4" & graph_values$t == "t+1"] <- "t+2" # Fix t-4
  graph_values$t[graph_values$name == "t-4" & graph_values$t == "t+0"] <- "t+1" # Fix t-4
  
  # Convert t to factor and fix order 
  graph_values$t <- as.factor(graph_values$t)
  levels(graph_values$t)
  graph_values$t_num <- as.numeric(graph_values$t)
  graph_values$name <- fct_rev(graph_values$name)
  
  # Plot Logged
  out <- graph_values %>% 
    filter(pooling != "Count") %>% 
    ggplot(aes(x = t_num, y = estimate)) +  
    scale_color_brewer(palette= "Dark2") + 
    facet_wrap(~ name, nrow=1) + 
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  width=0.075, position = position_dodge(width = 0.75)) +
    geom_crossbar(aes(ymin = lower_90, ymax = upper_90), 
                  width=0.05, position = position_dodge(width = 0.75), color = "black") +
    geom_point(aes(x = t_num, y = estimate),
               size = 3, 
               position = position_dodge(width = 0.75)) + 
    scale_shape_manual(values=c(1, 16)) + 
    theme_bw() +
    geom_hline(yintercept=0, linetype="dashed", color = "black") +
    scale_x_continuous(breaks=c(1,3,5),
                       labels=c("t-5", "t-3", "t-1")) + 
    xlab("Years relative to placebo withdrawal") + 
    ylab("Percent change\nin reforms (ATT)") +
    theme(legend.position= "none", 
          legend.title=element_blank(), 
          strip.text = element_text(size = 14),
          legend.text = element_text(size=14),
          plot.title = element_text(hjust = 0.5, size=16), 
          text = element_text(size = 10),
          axis.title = element_blank(), 
          axis.text = element_text(size = 14))
  return(out)
}

# Supp. Placebo Plots -----

# Direct Placebos
placebo_plot_fxn(ests_cbps_p1_direct, ests_cbps_p2_direct, ests_cbps_p3_direct,
                 ests_cbps_p4_direct, ests_cbps_p5_direct)
ggsave("figs/placebo_direct.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_rel_power, ests_cbps_p2_direct_rel_power, ests_cbps_p3_direct_rel_power,
                 ests_cbps_p4_direct_rel_power, ests_cbps_p5_direct_rel_power)
ggsave("figs/placebo_direct_power.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_not_rel_power, ests_cbps_p2_direct_not_rel_power, ests_cbps_p3_direct_not_rel_power,
                 ests_cbps_p4_direct_not_rel_power, ests_cbps_p5_direct_not_rel_power)
ggsave("figs/placebo_direct_not_power.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile1, ests_cbps_p2_direct_tertile1, ests_cbps_p3_direct_tertile1,
                 ests_cbps_p4_direct_tertile1, ests_cbps_p5_direct_tertile1)
ggsave("figs/placebo_direct_tertile1.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile2, ests_cbps_p2_direct_tertile2, ests_cbps_p3_direct_tertile2,
                 ests_cbps_p4_direct_tertile2, ests_cbps_p5_direct_tertile2)
ggsave("figs/placebo_direct_tertile2.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile3, ests_cbps_p2_direct_tertile3, ests_cbps_p3_direct_tertile3,
                 ests_cbps_p4_direct_tertile3, ests_cbps_p5_direct_tertile3)
ggsave("figs/placebo_direct_tertile3.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

# Indirect Placebos
placebo_plot_fxn(ests_cbps_p1_indirect, ests_cbps_p2_indirect, ests_cbps_p3_indirect,
                 ests_cbps_p4_indirect, ests_cbps_p5_indirect)
ggsave("figs/placebo_indirect.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_indirect_power, ests_cbps_p2_indirect_power, ests_cbps_p3_indirect_power,
                 ests_cbps_p4_indirect_power, ests_cbps_p5_indirect_power)
ggsave("figs/placebo_indirect_power.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_indirect_not_power, ests_cbps_p2_indirect_not_power, ests_cbps_p3_indirect_not_power,
                 ests_cbps_p4_indirect_not_power, ests_cbps_p5_indirect_not_power)
ggsave("figs/placebo_indirect_not_power.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_indirect_tertile1, ests_cbps_p2_indirect_tertile1, ests_cbps_p3_indirect_tertile1,
                 ests_cbps_p4_indirect_tertile1, ests_cbps_p5_indirect_tertile1)
ggsave("figs/placebo_indirect_tertile1.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_indirect_tertile2, ests_cbps_p2_indirect_tertile2, ests_cbps_p3_indirect_tertile2,
                 ests_cbps_p4_indirect_tertile2, ests_cbps_p5_indirect_tertile2)
ggsave("figs/placebo_indirect_tertile2.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_indirect_tertile3, ests_cbps_p2_indirect_tertile3, ests_cbps_p3_indirect_tertile3,
                 ests_cbps_p4_indirect_tertile3, ests_cbps_p5_indirect_tertile3)
ggsave("figs/placebo_indirect_tertile3.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

# Direct Restricted Placebos
placebo_plot_fxn(ests_cbps_p1_direct_restricted, ests_cbps_p2_direct_restricted, ests_cbps_p3_direct_restricted,
                 ests_cbps_p4_direct_restricted, ests_cbps_p5_direct_restricted)
ggsave("figs/placebo_direct_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_rel_power_restricted, ests_cbps_p2_direct_rel_power_restricted, ests_cbps_p3_direct_rel_power_restricted,
                 ests_cbps_p4_direct_rel_power_restricted, ests_cbps_p5_direct_rel_power_restricted)
ggsave("figs/placebo_direct_power_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_not_rel_power_restricted, ests_cbps_p2_direct_not_rel_power_restricted, ests_cbps_p3_direct_not_rel_power_restricted,
                 ests_cbps_p4_direct_not_rel_power_restricted, ests_cbps_p5_direct_not_rel_power_restricted)
ggsave("figs/placebo_direct_not_power_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile1_restricted, ests_cbps_p2_direct_tertile1_restricted, ests_cbps_p3_direct_tertile1_restricted,
                 ests_cbps_p4_direct_tertile1_restricted, ests_cbps_p5_direct_tertile1_restricted)
ggsave("figs/placebo_direct_tertile1_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile2_restricted, ests_cbps_p2_direct_tertile2_restricted, ests_cbps_p3_direct_tertile2_restricted,
                 ests_cbps_p4_direct_tertile2_restricted, ests_cbps_p5_direct_tertile2_restricted)
ggsave("figs/placebo_direct_tertile2_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")

placebo_plot_fxn(ests_cbps_p1_direct_tertile3_restricted, ests_cbps_p2_direct_tertile3_restricted, ests_cbps_p3_direct_tertile3_restricted,
                 ests_cbps_p4_direct_tertile3_restricted, ests_cbps_p5_direct_tertile3_restricted)
ggsave("figs/placebo_direct_tertile3_restricted.jpeg", plot = last_plot(), width = 6.5, height = 2, units = "in")


# Treatment and Outcome Distributions ----

# Treatment distributions
DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "direct_exit_initial", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_direct.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "direct_exit_initial_not_rel_power", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_direct_not_power.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "direct_exit_initial_rel_power", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_direct_power.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "indirect_exit_initial", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_indirect.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "indirect_exit_initial_not_power", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_indirect_not_power.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

DisplayTreatment(unit.id = "regime_number",
                 time.id = "year", legend.position = "none",
                 title = element_blank(),
                 xlab = "Year", ylab = "Regime",
                 treatment = "indirect_exit_initial_power", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/treatment_dist_indirect_power.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")

# Outcome distribution
analysis_dat <- analysis_dat |> mutate(reform_indicator = as.integer(reform_n > 0))
DisplayTreatment(unit.id = "regime_number",
                 title = element_blank(),
                 time.id = "year", legend.position = "none",
                 xlab = "Year", ylab = "Regime",
                 treatment = "reform_indicator", data = analysis_dat,
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE)
ggsave("figs/outcomes_dist.jpeg", plot = last_plot(), width = 4, height = 3, units = "in")


# Map of Exits -----

# Load data
world <- map_data("world")
states <- read.csv("data/preprocessing/states2016.csv")

# Clean and merge
world <- world %>% filter(region != "Antarctica")
map_countries <- unique(world$region)

# Check overlap
in_map <- subset(states$statenme, states$statenme %in% map_countries)
not_in_map <- subset(states$statenme, !(states$statenme %in% map_countries))

# Check overlap
not_in_cow <- subset(map_countries, !(map_countries %in% states$statenme))
in_cow <- subset(map_countries, (map_countries %in% states$statenme))

# Fix naming (Changes not_in_cow name to not_in_map name)
world$region[world$region == "USA"] <- "United States of America"
world$region[world$region == "UK"] <- "United Kingdom"
world$region[world$region == "Timor-Leste"] <- "East Timor"
world$region[world$region == "Republic of Congo"] <- "Congo"
world$region[world$region == "Serbia"] <- "Yugoslavia"
world$region[world$region == "Micronesia"] <- "Federated States of Micronesia"

# Merge
world$statenme <- world$region
world <- left_join(world, states, by = "statenme", relationship = "many-to-many")

# Clean treaty ratification data 
map_dat <- unilateral_exits %>% group_by(part_ccode) %>% summarize(
  wds = n())
map_dat$year <- 2016
colnames(map_dat)[1] <- "ccode"

# Merge and fill in no withdrawal states with 0 
world <- left_join(world, map_dat, by = "ccode")
world[c("wds")][is.na(world[c("wds")])] <- 0

# Reformat for plotting
vars <- c("wds") # "rev_troika_in_current_session", "treatment_sessions")
long_world <- pivot_longer(world, cols = all_of(vars))
long_world$label <- long_world$name
long_world$label[long_world$name == "wds"] <- "Initial Unilateral Exits"
Xs <- c("Initial Unilateral Exits")
long_world$label <- factor(long_world$label, levels = Xs)
long_world_backup <- long_world

# Bin values for legend Updated
long_world$value[long_world$value > 0 & long_world$value <= 1] <- 1 
long_world$value[long_world$value > 1 & long_world$value <= 5] <- 2 
long_world$value[long_world$value > 5] <- 3
long_world$value <- factor(long_world$value)

# Make map
subset(long_world, label == "Initial Unilateral Exits") %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = value)) +
  geom_polygon(color = "black", linewidth = 0.1, alpha = 1.0) + 
  scale_fill_manual(name="Initial Unilateral Exits: ",
                    values=c("lightgrey", "lightblue1", "lightblue3", "lightblue4", "lightblue4"),
                    labels=c("0", "1", "2-5", ">5")) +
  theme_bw() +
  theme(legend.position= "bottom", 
        legend.title=element_text(size=16, color = "black"), #"#002E63"), 
        legend.text = element_text(size=16, color = "black"),
        panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.border = element_blank()) 
ggsave("figs/initial_unilateral_exit_map.jpeg", # Save plot
       plot =  last_plot(), width = 8.5, height = 5, units = "in") 

# Secret Alliances ----
# Import and clean data
atop <- read.csv("data/preprocessing/atop5_0a.csv", header=T) # Import ATOP Data
atop_dat <- subset(atop, atop$pubsecr==2) # Subset to secret alliances

# Assign beginning year to one treaty without end year
atop_dat$endyr[atop_dat$endyr == 0] <- atop_dat$begyr[atop_dat$endyr == 0]
sec_treaty_years <- numeric() # Loop years into vector
for(i in 1:nrow(atop_dat)){sec_treaty_years <- c(sec_treaty_years, atop_dat$begyr[i]:atop_dat$endyr[i])}

# Format data for graphing
year_dat <- as.data.frame(table(sec_treaty_years)) # Coerce to dataframe 
colnames(year_dat) <- c("Year", "Count") # Assign column names 
year_dat$Year <- as.numeric(as.character(year_dat$Year)) # Fix formatting

# Fix that some years are not represented
vals <- as.data.frame(matrix(NA, nrow = 204, ncol = 2))
colnames(vals) <- c("Year", "Count_fill")
vals$Year <- 1815:2018
vals$Count_fill <- 0
year_dat <- merge(vals, year_dat, by = "Year", all = T)
year_dat$Count[is.na(year_dat$Count)] <- year_dat$Count_fill[is.na(year_dat$Count)] # Fill in

ggplot(year_dat, aes(x=Year,y=Count)) +
  geom_line() +
  theme_bw() + labs(x = "Year", y = "Number of Secret Alliances") +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_vline(xintercept = 1918, color = "red", linetype = "dashed") +
  annotate("text",x=c(1920),y=14,label=c("1918"),hjust=0, color = "red") +
  scale_x_continuous(minor_breaks = seq(1800, 2020, 25), breaks = seq(1800, 2000, 25))
ggsave("figs/secret_alliances_per_year.jpeg", plot = last_plot(), width = 8, height = 4, units = "in")


# Substantive Results ----

ests_cbps_direct$estimates # Direct effect 
1/mean(analysis_dat$reform_n, na.rm = T) # Average number of reforms
1/mean(analysis_dat$reform_n[analysis_dat$direct_exit_initial_placebo_1 == 1], na.rm = T) # Average number of reforms pre-withdrawal
mean(analysis_dat$reform_n[analysis_dat$direct_exit_initial_placebo_1 == 1], na.rm = T) * length(ests_cbps_direct$estimates) # Expected number of reforms without withdrawal 
sum(ests_cbps_direct$estimates) # Additional reforms after exit 
sum(ests_cbps_direct$estimates) + mean(analysis_dat$reform_n[analysis_dat$direct_exit_initial_placebo_1 == 1], na.rm = T) * length(ests_cbps_direct$estimates) # Expected number of reforms with withdrawal 
sum(ests_cbps_direct$estimates) * sum(analysis_dat$direct_exit_initial) # Total effect across study period
ests_cbps_direct_rel_power$estimates # Point estimates of major power withdrawal 
ests_cbps_direct_rel_power$estimates > ests_cbps_direct$estimates # All major power estimates are greater than main estimates 
ests_cbps_direct_not_rel_power$estimates > ests_cbps_direct$estimates # All minor power estimates are less than main estimates 
sum(ests_cbps_direct_rel_power$estimates) * sum(analysis_dat$direct_exit_initial_rel_power) # Total effect across study period (major powers)
sum(ests_cbps_direct_not_rel_power$estimates) * sum(analysis_dat$direct_exit_initial_not_rel_power) # Total effect across study period (minor powers)

ests_cbps_indirect$estimates # Indirect effect 
mean(analysis_dat$indirect_exit_initial) # Proprotion of units treated indirectly
mean(analysis_dat$direct_exit_initial) # Proprotion of units treated directly
sum(ests_cbps_indirect$estimates) * sum(analysis_dat$indirect_exit_initial) # Substantive effect 
(sum(ests_cbps_indirect$estimates) * sum(analysis_dat$indirect_exit_initial)) + (sum(ests_cbps_direct$estimates) * sum(analysis_dat$direct_exit_initial)) # Net effect 
ests_cbps_indirect_power$estimates < ests_cbps_indirect$estimates # Comparison of indirect effect power estimates to pooled indirect effect estimates
ests_cbps_indirect_not_power$estimates < ests_cbps_indirect$estimates
sum(ests_cbps_indirect_power$estimates) * sum(analysis_dat$indirect_exit_initial_power) # Substantive effects 
sum(ests_cbps_indirect_not_power$estimates) * sum(analysis_dat$indirect_exit_initial_not_power)
apply(ests_cbps_direct$bootstrapped.estimates, 2, high_ci_lower_bound_fxn) + ests_cbps_indirect$estimates # Note positive lower bound of 95%-CI of direct effect after subtracting indirect effect estimates

# Increase in percentage of international laws had no withdrawals occured 1945-2018
regimes <- read.csv("data/treaty_data/regime_ids.csv") # Import treaty-regime key
n_regimes <- sum(regimes$amendment_code == 1 | regimes$core_treaty == 1) # Total reforms and parent agreements
total_direct_effect <- sum(ests_cbps_direct$estimates) * sum(analysis_dat$direct_exit_initial) # Total effects
total_indirect_effect <- sum(ests_cbps_indirect$estimates) * sum(analysis_dat$indirect_exit_initial) 
combined_effect <- total_direct_effect + total_indirect_effect # Combined effect
(n_regimes + (-1 * combined_effect))/n_regimes # Remove combined effect from observed total, compute percentage change

# Summary Stats ----

# Number regimes
length(unique(analysis_dat$regime_number))

# Make list of relevant variables
vars_desc <- c(
  
  # Outcomes
  "reform_n",
  
  # Treatment
  "direct_exit_initial",
  "direct_exit_initial_rel_power",
  "direct_exit_initial_not_rel_power",
  "indirect_exit_initial",
  "indirect_exit_initial_power",
  "indirect_exit_initial_not_power",
  
  # Network Variables             
  "overlap", 
  "overlap_rel_power",
  "overlap_not_rel_power",
  "overlap_power",
  "overlap_not_power",
  
  # Covariates
  "n_members_log", 
  "last_reform", 
  "cumulative_agreements",
  "pref_mean", 
  "pref_sd",
  #"human_rights_count", 
  "human_rights", 
  #"environment_count",
  "environment",
  #"security_count", 
  "security", 
  #"economics_count",
  "economics",
  
  # Index
  "regime_number_int",
  "year"
  
)

# Make list of variable names
var_names <- c(
  
  # Outcomes
  "Reforms (Count)", # "reform_n",
  
  # Treatment
  "Direct Exit (All)", #  "direct_exit_initial",
  "Direct Exit (Major Power)", #  "direct_exit_initial_power",
  "Direct Exit (Minor Power)", #  "direct_exit_initial_not_power",
  "Indirect Exit (All)", #  "indirect_exit_initial",
  "Indirect Exit (Major Power)", #  "indirect_exit_initial_power",
  "Indirect Exit (Minor Power)", #  "indirect_exit_initial_not_power",
  
  # Network Variables             
  "Overlap (All)", # "overlap", 
  "Overlap (Major Power, Direct)", # "overlap_rel_power",
  "Overlap (Minor Power, Direct)", # "overlap_not_rel_power",
  "Overlap (Major Power, Indirect)", # "overlap_power",
  "Overlap (Minor Power, Indirect)", # "overlap_not_power",
  
  # Covariates
  "Members (Log)",  # "n_members_log", 
  "Years Since Last Reform",  #  "last_reform", 
  "Cumulative Reforms",  # "cumulative_agreements",
  "UNGA Ideal Point (Mean)",  # "pref_mean", 
  "UNGA Ideal Point (St. Dev.)",  # "pref_sd"
  "Human Rights (Count)",  #  "human_rights", 
  "Environment (Count)",  #  "environment", 
  "Security (Count)",  #  "security", 
  "Economics (Count)",  #  "economics",
  
  # Index
  "Regime ID",  #  "regime_number_int",
  "Year"  #  "year"
)

# Subset relevant variables for full sample
temp_data <- analysis_dat[,vars_desc]
footnote <- c("To be replaced.") # A placeholder
desc_stat_full <- stargazer(as.data.frame(temp_data), digits = 2,
                            #type = "text",
                            float = FALSE,
                            notes = footnote,
                            covariate.labels = var_names)

# Fix note formatting
note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{\\textit{Note:} Summary statistics computed using all observations in full dataset; data used in analyses differs because of matching and weighting. }} \\\\"
desc_stat_full[grepl("To be replaced.", desc_stat_full)] <- note.latex
cat(desc_stat_full, sep = "\n")

# Add resizebox and float to latex output
desc_stat_full <- c("\\begin{table}[!htbp] \\centering",
                    "\\caption{Descriptive Statistics}",
                    "\\label{desc_stat_full}",
                    desc_stat_full,
                    "\\end{table}")

# Save
file_connection <-file("figs/desc_stat.txt")
writeLines(desc_stat_full, file_connection)
close(file_connection)


# Matched Set Examples -----

# Direct Effect Analysis
direct_weights <- numeric()
for(i in 1:length(matches_cbps_direct$att)){
  direct_weights <- c(direct_weights, unname(attr(matches_cbps_direct$att[[i]], "weights")))
}
summary(matches_cbps_direct$att)$set.size.summary # Summary of matched set size
summary(direct_weights) # Summary of weights
hist(direct_weights, breaks = 1000) # Histograms of weights

# Indirect Effect Analysis
indirect_weights <- numeric()
for(i in 1:length(matches_cbps_indirect$att)){
  indirect_weights <- c(indirect_weights, unname(attr(matches_cbps_indirect$att[[i]], "weights")))
}
summary(matches_cbps_indirect$att)$set.size.summary  # Summary of matched set size
summary(indirect_weights) # Summary of weights
hist(indirect_weights, breaks = 1000) # Histograms of weights

# Top Matches in Direct Effect Analysis
rm(table_data)
check <- summary(matches_cbps_direct$att)$overview
set.seed(123) # 123 is good
weighting_examples <- sample(1:nrow(check), 10, replace = F)
for(i in 1:length(weighting_examples)){
  j <- weighting_examples[i] # Extract row number
  # Get to match regime number and weight
  top_matched_regime <- attr(matches_cbps_direct$att[[j]], "weights")[as.numeric(which.max(attr(matches_cbps_direct$att[[j]], "weights")))]
  table_data_temp <- tibble( # Assign values to table for printing
    "Treated Regime" = analysis_dat$title_of_max_rat[
      analysis_dat$regime_number_int == check$regime_number_int[j] & 
        analysis_dat$year == check$year[j]],
    "Year" = check$year[j],
    "Matched Set Size" = check$matched.set.size[j],
    "Most Weighted Control Regime" = analysis_dat$title_of_max_rat[
      analysis_dat$regime_number_int == as.numeric(names(top_matched_regime)) & 
        analysis_dat$year == check$year[j]],
    "Weight" = unname(top_matched_regime)
  )
  if(exists("table_data")){
    table_data <- rbind(table_data, table_data_temp)
  }else{
    table_data <- table_data_temp
  }
}

# Make Table
table_out <- print(
  xtable( # Table contents
    table_data[1:5,],
    align=c(
      "|L{0\\textwidth}|",
      "|L{0.35\\textwidth}|",
      "L{0.05\\textwidth}|",
      "L{0.075\\textwidth}|",
      "L{0.35\\textwidth}|",
      "L{0.075\\textwidth}|"),
    caption = "\\emph{Matched Regimes with Most Weight, Direct Effect}. Examples of regimes assigned most weight to randomly selected cases of regimes directly treated by withdrawal."),
  booktabs = TRUE, # Format table
  size="\\fontsize{10pt}{10pt}\\selectfont",
  include.rownames = F,
  hline.after = -1:5# rep(c(1:nrow(table_data)),2)
)

# Save
label <- "top_matches_direct" 
file_path_name <- str_c(c("figs/", label, ".txt"), collapse =  "")
file_connection <- file(file_path_name) # Save as txt in figs folder
writeLines(table_out, file_connection)
close(file_connection)

#' Note: regime names in top_matches_direct.txt are of the treaty in each regime with the most ratifications. 
#' To make the table more comprehensible I manually edit regime names so that common names are used. 
#' For example, changing "(a) Articles of Agreement of the International Monetary Fund" to "International Monitary Fund"
#' This is then saved as top_matches_direct_revised_names.txt 

# Top Matches in Indirect Effect Analysis
rm(table_data)
check <- summary(matches_cbps_indirect$att)$overview
weighting_examples <- sample(1:nrow(check), 10, replace = F)
for(i in 1:length(weighting_examples)){
  j <- weighting_examples[i] # Extract row number
  # Get to match regime number and weight
  top_matched_regime <- attr(matches_cbps_indirect$att[[j]], "weights")[as.numeric(which.max(attr(matches_cbps_indirect$att[[j]], "weights")))]
  table_data_temp <- tibble( # Assign values to table for printing
    "Treated Regime" = analysis_dat$title_of_max_rat[
      analysis_dat$regime_number_int == check$regime_number_int[j] & 
        analysis_dat$year == check$year[j]],
    "Year" = check$year[j],
    "Matched Set Size" = check$matched.set.size[j],
    "Most Weighted Control Regime" = analysis_dat$title_of_max_rat[
      analysis_dat$regime_number_int == as.numeric(names(top_matched_regime)) & 
        analysis_dat$year == check$year[j]],
    "Weight" = unname(top_matched_regime)
  )
  if(exists("table_data")){
    table_data <- rbind(table_data, table_data_temp)
  }else{
    table_data <- table_data_temp
  }
}

# Make Table
table_out <- print(
  xtable( # Table contents
    table_data[1:5,],
    align=c(
      "|L{0\\textwidth}|",
      "|L{0.35\\textwidth}|",
      "L{0.05\\textwidth}|",
      "L{0.075\\textwidth}|",
      "L{0.35\\textwidth}|",
      "L{0.075\\textwidth}|"),
    caption = "\\emph{Matched Regimes with Most Weight, Indirect Effect}. Examples of regimes assigned most weight to randomly selected cases of regimes indirectly treated by withdrawal."),
  booktabs = TRUE, # Format table
  size="\\fontsize{10pt}{10pt}\\selectfont", 
  include.rownames = F,
  hline.after = -1:5# rep(c(1:nrow(table_data)),2)
)

# Save
label <- "top_matches_indirect" 
file_path_name <- str_c(c("figs/", label, ".txt"), collapse =  "")
file_connection <- file(file_path_name) # Save as txt in figs folder
writeLines(table_out, file_connection)
close(file_connection)

#' Note: regime names in top_matches_indirect.txt are of the treaty in each regime with the most ratifications. 
#' To make the table more comprehensible I manually edit regime names so that common names are used. 
#' For example, changing "(a) Articles of Agreement of the International Monetary Fund" to "International Monitary Fund"
#' This is then saved as top_matches_indirect_revised_names.txt 
